home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / msdos / raytrace / pov / source / image.c < prev    next >
C/C++ Source or Header  |  1993-07-28  |  22KB  |  762 lines
  1. /****************************************************************************
  2. *                image.c
  3. *
  4. *  This module implements the mapped textures including image map, bump map
  5. *  and material map. 
  6. *
  7. *  from Persistence of Vision Raytracer
  8. *  Copyright 1993 Persistence of Vision Team
  9. *---------------------------------------------------------------------------
  10. *  NOTICE: This source code file is provided so that users may experiment
  11. *  with enhancements to POV-Ray and to port the software to platforms other 
  12. *  than those supported by the POV-Ray Team.  There are strict rules under
  13. *  which you are permitted to use this file.  The rules are in the file
  14. *  named POVLEGAL.DOC which should be distributed with this file. If 
  15. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  16. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  17. *  Forum.  The latest version of POV-Ray may be found there as well.
  18. *
  19. * This program is based on the popular DKB raytracer version 2.12.
  20. * DKBTrace was originally written by David K. Buck.
  21. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  22. *
  23. *****************************************************************************/
  24.  
  25. #include "frame.h"
  26. #include "vector.h"
  27. #include "povproto.h"
  28. #include "texture.h"
  29.  
  30. static int cylindrical_image_map PARAMS((DBL x, DBL y, DBL z, IMAGE *Image, DBL *u, DBL *v));
  31. static int torus_image_map PARAMS((DBL x, DBL y, DBL z, IMAGE *Image, DBL *u, DBL *v));
  32. static int spherical_image_map PARAMS((DBL x, DBL y, DBL z, IMAGE *Image, DBL *u, DBL *v));
  33. static int planar_image_map PARAMS((DBL x, DBL y, DBL z, IMAGE *Image, DBL *u, DBL *v));
  34. static void no_interpolation PARAMS((IMAGE *Image,DBL xcoor,DBL ycoor,COLOUR *colour,int *index));
  35. static DBL bilinear PARAMS((DBL *corners,DBL x,DBL y));
  36. static DBL norm_dist PARAMS((DBL *corners,DBL x,DBL y));
  37. static void Interp PARAMS((IMAGE *Image,DBL xcoor,DBL ycoor,COLOUR *colour,int *index));
  38. static void image_colour_at PARAMS((IMAGE *Image,DBL xcoor,DBL ycoor,COLOUR *colour,int *index));
  39.  
  40. /*
  41.    2-D to 3-D Procedural Texture Mapping of a Bitmapped Image onto an Object:
  42.            
  43. A. Simplistic (planar) method of image projection devised by DKB and AAC:
  44.  
  45.    1. Transform texture in 3-D space if requested.
  46.    2. Determine local object 2-d coords from 3-d coords by <X Y Z> triple.
  47.    3. Return pixel color value at that position on the 2-d plane of "Image".
  48.    3. Map colour value in Image [0..255] to a more normal colour range [0..1].
  49.  
  50. B. Specialized shape projection variations by Alexander Enzmann:
  51.  
  52.    1. Cylindrical mapping
  53.    2. Spherical mapping
  54.    3. Torus mapping
  55. */
  56.  
  57. void image_map (x, y, z, Pigment, colour)
  58. DBL x, y, z;
  59. PIGMENT *Pigment;
  60. COLOUR *colour;
  61.   {
  62.   /* determine local object 2-d coords from 3-d coords */
  63.   /* "unwrap" object 2-d coord onto flat 2-d plane */
  64.   /* return pixel color value at that posn on 2-d plane */
  65.  
  66.   DBL xcoor = 0.0, ycoor = 0.0    ;
  67.   int reg_number;
  68.  
  69.   if(map (x, y, z, ((TPATTERN *)Pigment), Pigment->Image,&xcoor, &ycoor))
  70.     {
  71.     Make_ColourA (colour,1.0,1.0,1.0,1.0);
  72.     return;
  73.     }
  74.   image_colour_at(Pigment->Image,xcoor, ycoor, colour, ®_number);
  75.   }
  76.  
  77. /* Very different stuff than the other routines here. This routine takes  */
  78. /* an intersection point and a texture and returns a new texture based on */
  79. /* the index/color of that point in an image/materials map. CdW 7/91      */
  80.  
  81. TEXTURE *material_map(IPoint,Texture)
  82. VECTOR *IPoint;
  83. MATERIAL *Texture;
  84.   {
  85.   register DBL x, y, z;
  86.   DBL xcoor = 0.0, ycoor = 0.0;
  87.   int reg_number = 0;
  88.   COLOUR colour;
  89.   int Material_Number=0;
  90.   TEXTURE *Temp_Tex;
  91.   int numtex;
  92.   VECTOR TPoint;
  93.   TPATTERN TPattern;
  94.  
  95.   INIT_TPATTERN_FIELDS((&TPattern),0);
  96.   Make_ColourA(&colour,0.0,0.0,0.0,0.0);
  97.  
  98.   if (Texture->Trans != NULL) 
  99.     MInvTransPoint (&TPoint, IPoint, Texture->Trans);
  100.   else 
  101.     TPoint = *IPoint;
  102.  
  103.   x = TPoint.x;
  104.   y = TPoint.y;
  105.   z = TPoint.z;
  106.  
  107.   /* now we have transformed x, y, z we use image mapping routine */
  108.   /* to determine texture index */
  109.   if(map (x, y, z, &TPattern, Texture->Image,&xcoor, &ycoor))
  110.     Material_Number=0;
  111.     else{
  112.     image_colour_at(Texture->Image,xcoor, ycoor, &colour, ®_number);
  113.  
  114.     if (Texture->Image->Colour_Map == NULL) 
  115.       Material_Number = (int)colour.Red*255;
  116.     else 
  117.       Material_Number = reg_number;
  118.     } 
  119.  
  120.   if(Material_Number > Texture->Num_Of_Mats)
  121.     Material_Number %= Texture->Num_Of_Mats;
  122.   for(numtex=0, Temp_Tex = Texture->Materials;
  123.   (Temp_Tex->Next_Material != NULL) && (numtex<Material_Number); 
  124.   Temp_Tex = Temp_Tex->Next_Material, numtex++)
  125.     ; /* do nothing */
  126.  
  127.   return(Temp_Tex);  
  128.   }
  129.  
  130. TEXTURE *tiles_texture(IPoint,Texture)
  131. VECTOR *IPoint;
  132. TILES *Texture;
  133.   {
  134.   int Block;
  135.   VECTOR TPoint;
  136.  
  137.   if (Texture->Trans != NULL) 
  138.     MInvTransPoint (&TPoint, IPoint, Texture->Trans);
  139.   else 
  140.     TPoint = *IPoint;
  141.  
  142.   Block = (int)(FLOOR(TPoint.x+Small_Tolerance)
  143.     +FLOOR(TPoint.y+Small_Tolerance)
  144.     +FLOOR(TPoint.z+Small_Tolerance));
  145.  
  146.   if (Block & 1)
  147.     return(Texture->Tile1);
  148.   else
  149.     return(Texture->Tile2);
  150.  
  151.   } 
  152.  
  153.  
  154. void bump_map (x, y, z, Tnormal, normal)
  155. DBL x, y, z;
  156. TNORMAL *Tnormal;
  157. VECTOR *normal;
  158.   {
  159.   DBL xcoor = 0.0, ycoor = 0.0;
  160.   int index,index2,index3;
  161.   COLOUR colour, colour2, colour3;
  162.   VECTOR p1,p2,p3;
  163.   VECTOR bump_normal;
  164.   VECTOR xprime, yprime, zprime, Temp;
  165.   DBL Length;
  166.   DBL Amount = Tnormal->Amount;
  167.   IMAGE *Image = Tnormal->Image;
  168.   Make_ColourA (&colour,0.0,0.0,0.0,0.0);
  169.   Make_ColourA (&colour2,0.0,0.0,0.0,0.0);
  170.   Make_ColourA (&colour3,0.0,0.0,0.0,0.0);
  171.  
  172.   /* going to have to change this */
  173.   /* need to know if bump point is off of image for all 3 points */
  174.  
  175.   if(map (x, y, z, (TPATTERN *)Tnormal, Image,&xcoor, &ycoor))
  176.     return;
  177.   else
  178.     image_colour_at(Image,xcoor, ycoor, &colour, &index);
  179.  
  180.   xcoor--;
  181.   ycoor++;
  182.   if (xcoor < 0.0)
  183.     xcoor += (DBL)Image->iwidth;
  184.   else if (xcoor >= Image->iwidth)
  185.     xcoor -= (DBL)Image->iwidth;
  186.   if (ycoor < 0.0)
  187.     ycoor += (DBL)Image->iheight;
  188.   else if (ycoor >=(DBL) Image->iheight)
  189.     ycoor -= (DBL)Image->iheight;
  190.   image_colour_at(Image,xcoor, ycoor, &colour2, &index2);
  191.  
  192.   xcoor +=2.0;
  193.   if (xcoor < 0.0)
  194.     xcoor += (DBL)Image->iwidth;
  195.   else if (xcoor >= Image->iwidth)
  196.     xcoor -= (DBL)Image->iwidth;
  197.  
  198.   image_colour_at(Image,xcoor, ycoor, &colour3, &index3);
  199.  
  200.   if (Options & DEBUGGING)
  201.     printf ("Bump Map %g %g %g xcoor %g ycoor %g\n", x, y, z, xcoor, ycoor);
  202.  
  203.   if (Image->Colour_Map == NULL || 
  204.     Image->Use_Colour_Flag ) 
  205.     {
  206.     p1.x = 0; 
  207.     p1.y =
  208.     Amount*(0.229*colour.Red+0.587*colour.Green+0.114*colour.Blue);
  209.     p1.z = 0;
  210.     p2.x = 0; 
  211.     p2.y =
  212.     Amount*(0.229*colour2.Red+0.587*colour2.Green+0.114*colour2.Blue);
  213.     p2.z = 1;
  214.     p3.x = 1; 
  215.     p3.y =
  216.     Amount*(0.229*colour3.Red+0.587*colour3.Green+0.114*colour3.Blue);
  217.     p3.z = 1;
  218.     }
  219.   else 
  220.     {
  221.     p1.x = 0; 
  222.     p1.y = Amount * index;
  223.     p1.z = 0;
  224.     p2.x = 0; 
  225.     p2.y = Amount * index2;
  226.     p2.z = 1;
  227.     p3.x = 1; 
  228.     p3.y = Amount * index3;
  229.     p3.z = 1;
  230.     }
  231.   /* we have points 1,2,3 for a triangle now we need the surface normal for it
  232. */  
  233.   VSub (xprime, p1, p2);
  234.   VSub (yprime, p3, p2);
  235.   VCross (bump_normal, yprime, xprime);
  236.   VNormalize(bump_normal, bump_normal);
  237.  
  238.     Make_Vector(&yprime,normal->x,normal->y,normal->z);
  239.   Make_Vector(&Temp,0.0,1.0,0.0);
  240.   VCross (xprime,yprime,Temp);
  241.   VLength(Length,xprime);
  242.   if(Length < 1.0e-9)
  243.     {
  244.     if(fabs(normal->y - 1.0)<Small_Tolerance)
  245.       {
  246.       Make_Vector(&yprime,0.0,1.0,0.0);
  247.       Make_Vector(&xprime,1.0,0.0,0.0);
  248.       Length = 1.0;
  249.       }
  250.     else
  251.       {
  252.       Make_Vector(&yprime,0.0,-1.0,0.0);
  253.       Make_Vector(&xprime,1.0,0.0,0.0);
  254.       Length = 1.0;
  255.       }
  256.     }
  257.   VScale(xprime,xprime,1.0/Length);
  258.   VCross (zprime,xprime,yprime);
  259.   VNormalize(zprime,zprime);
  260.   VScale(xprime,xprime,bump_normal.x);
  261.   VScale(yprime,yprime,bump_normal.y);
  262.   VScale(zprime,zprime,bump_normal.z);
  263.   VAdd(Temp,xprime,yprime);
  264.   VAdd(*normal,Temp,zprime);
  265.   VNormalize(*normal,*normal);
  266.  
  267.   return;
  268.   }
  269.  
  270.  
  271. static void image_colour_at(Image, xcoor, ycoor, colour, index)
  272. IMAGE *Image;     
  273. DBL xcoor, ycoor;
  274. COLOUR *colour;
  275. int *index;
  276.   {  
  277.  
  278.   switch(Image->Interpolation_Type)
  279.   {
  280.   case NO_INTERPOLATION:
  281.     no_interpolation(Image,xcoor,ycoor,colour,index);
  282.     break;
  283.   default: 
  284.     Interp(Image,xcoor, ycoor, colour, index);
  285.     break;
  286.   }   
  287.   }
  288.  
  289.  
  290.  
  291. /* Map a point (x, y, z) on a cylinder of radius 1, height 1, that has its
  292.    axis of symmetry along the y-axis to the square [0,1]x[0,1]. */
  293. static int cylindrical_image_map (x, y, z, Image, u, v)
  294. DBL x, y, z;
  295. IMAGE *Image;
  296. DBL *u, *v;
  297.   {
  298.   DBL len, theta;
  299.  
  300.   if ((Image->Once_Flag) && ((y < 0.0) || (y > 1.0)))
  301.     return 0;
  302.   *v = fmod (y * Image->height, Image->height);
  303.  
  304.   /* Make sure this vector is on the unit sphere. */
  305.   len = sqrt(x*x + y*y + z*z);
  306.   if (len == 0.0)
  307.     return 0;
  308.   else 
  309.     {
  310.     x /= len;
  311.     z /= len;
  312.     }   
  313.   /* Determine its angle from the point (1, 0, 0) in the x-z plane. */
  314.     len = sqrt(x*x + z*z);
  315.   if (len == 0.0)
  316.     return 0;
  317.   else 
  318.     {
  319.     if (z == 0.0)
  320.       if (x > 0)
  321.         theta = 0.0;
  322.       else
  323.         theta = M_PI;
  324.     else 
  325.       {
  326.       theta = acos(x / len);
  327.       if (z < 0.0) theta = 2.0 * M_PI - theta;
  328.       }
  329.     theta /= 2.0 * M_PI; /* This will be from 0 to 1 */
  330.     }
  331.   *u = (theta * Image->width);
  332.   return 1;
  333.   }
  334.  
  335. /* Map a point (x, y, z) on a torus  to a 2-d image. */
  336. static int torus_image_map (x, y, z, Image, u, v)
  337. DBL x, y, z;
  338. IMAGE *Image;
  339. DBL *u, *v;
  340.   {
  341.   DBL len, phi, theta;
  342.   DBL r0;
  343.  
  344.   r0 = Image->Gradient.x;
  345.  
  346.   /* Determine its angle from the x-axis. */
  347.   len = sqrt(x*x + z*z);
  348.   if (len == 0.0)
  349.     return 0;
  350.   else 
  351.     {
  352.     if (z == 0.0)
  353.       if (x > 0)
  354.         theta = 0.0;
  355.       else
  356.         theta = M_PI;
  357.     else 
  358.       {
  359.       theta = acos(x / len);
  360.       if (z < 0.0) theta = 2.0 * M_PI - theta;
  361.       }
  362.     }
  363.   theta = 0.0 - theta;
  364.  
  365.   /* Now rotate about the y-axis to get the point (x, y, z) into the x-y plane. */
  366.   x = len - r0;
  367.   len = sqrt(x * x + y * y);
  368.   phi = acos(-x / len);
  369.   if (y > 0.0) phi = 2.0 * M_PI - phi;
  370.  
  371.   /* Determine the parametric coordinates. */
  372.   theta /= 2.0 * M_PI;
  373.   phi /= 2.0 * M_PI;
  374.   *u = (-theta * Image->width);
  375.   *v = (phi * Image->height);
  376.   return 1;
  377.   }
  378.  
  379. /* Map a point (x, y, z) on a sphere of radius 1 to a 2-d image. (Or is it the
  380.    other way around?) */
  381. static int spherical_image_map (x, y, z, Image, u, v)
  382. DBL x, y, z;
  383. IMAGE *Image;
  384. DBL *u, *v;
  385.   {
  386.   DBL len, phi, theta;
  387.  
  388.   /* Make sure this vector is on the unit sphere. */
  389.   len = sqrt(x*x + y*y + z*z);
  390.   if (len == 0.0)
  391.     return 0;
  392.   else 
  393.     {
  394.     x /= len;
  395.     y /= len;
  396.     z /= len;
  397.     }
  398.   /* Determine its angle from the x-z plane. */
  399.     phi = 0.5 + asin(y) / M_PI; /* This will be from 0 to 1 */
  400.  
  401.   /* Determine its angle from the point (1, 0, 0) in the x-z plane. */
  402.   len = sqrt(x*x + z*z);
  403.   if (len == 0.0) 
  404.     {
  405.     /* This point is at one of the poles. Any value of xcoord will be ok...*/
  406.     theta = 0;
  407.     }
  408.   else 
  409.     {
  410.     if (z == 0.0)
  411.       if (x > 0)
  412.         theta = 0.0;
  413.       else
  414.         theta = M_PI;
  415.     else 
  416.       {
  417.       theta = acos(x / len);
  418.       if (z < 0.0) theta = 2.0 * M_PI - theta;
  419.       }
  420.     theta /= 2.0 * M_PI; /* This will be from 0 to 1 */
  421.     }
  422.   *u = (theta * Image->width);
  423.   *v = (phi * Image->height);
  424.   return 1;
  425.   }
  426.  
  427. /*
  428.    2-D to 3-D Procedural Texture Mapping of a Bitmapped Image onto an Object:
  429.        
  430.    Simplistic planar method of object image projection devised by DKB and AAC.
  431.  
  432.    1. Transform texture in 3-D space if requested.
  433.    2. Determine local object 2-d coords from 3-d coords by <X Y Z> triple.
  434.    3. Return pixel color value at that position on the 2-d plane of "Image".
  435.    3. Map colour value in Image [0..255] to a more normal colour range [0..1].
  436. */
  437.  
  438. /* Return 0 if there is no color at this point (i.e. invisible), return 1
  439.    if a good mapping is found. */
  440. static int planar_image_map(x, y, z, Image, u, v)
  441. DBL x, y, z;
  442. IMAGE *Image;
  443. DBL *u, *v;
  444.   {  
  445.   if (Image -> Gradient.x != 0.0) 
  446.     {
  447.     if ((Image->Once_Flag) &&
  448.       ((x < 0.0) || (x > 1.0)))
  449.       return 0;
  450.     if (Image -> Gradient.x > 0)
  451.       *u = fmod (x * Image->width, Image->width);
  452.     else
  453.       *v = fmod (x * Image->height, Image->height);
  454.     }
  455.   if (Image -> Gradient.y != 0.0) 
  456.     {
  457.     if ((Image->Once_Flag) &&
  458.       ((y < 0.0) || (y > 1.0)))
  459.       return 0;
  460.     if (Image -> Gradient.y > 0)
  461.       *u = fmod (y * Image->width, Image->width);
  462.     else
  463.       *v = fmod (y * Image->height, Image->height);
  464.     }
  465.   if (Image -> Gradient.z != 0.0) 
  466.     {
  467.     if ((Image->Once_Flag) &&
  468.       ((z < 0.0) || (z > 1.0)))
  469.       return 0;
  470.     if (Image -> Gradient.z > 0)
  471.       *u = fmod (z * Image->width, Image->width);
  472.     else
  473.       *v = fmod (z * Image->height, Image->height);
  474.     }
  475.   return 1;
  476.   }
  477. /* Map returns 1 if no color found (invisible) or 0 if color found */ 
  478. int map (x, y, z, TPattern, Image, xcoor, ycoor)
  479. DBL x, y, z;
  480. TPATTERN *TPattern;
  481. IMAGE *Image;
  482. DBL *xcoor, *ycoor;
  483.   {
  484.   /* determine local object 2-d coords from 3-d coords */
  485.   /* "unwrap" object 2-d coord onto flat 2-d plane */
  486.   /* return pixel color value at that posn on 2-d plane */
  487.  
  488.   /* This causes problems so let's do without it for this release
  489.    if (TPattern->Turbulence != 0.0)
  490.      {
  491.       DTurbulence (&TextureTurbulence, x, y, z,TPattern->Octaves);
  492.       x += TextureTurbulence.x * TPattern->Turbulence;
  493.       y += TextureTurbulence.y * TPattern->Turbulence;
  494.       z += TextureTurbulence.z * TPattern->Turbulence;
  495.      }
  496.    */
  497.  
  498.   /* Now determine which mapper to use. */
  499.   switch (Image->Map_Type) 
  500.   {
  501.   case PLANAR_MAP:
  502.     if (!planar_image_map(x, y, z, Image, xcoor, ycoor)) 
  503.       return(1);
  504.     break;
  505.   case SPHERICAL_MAP:
  506.     if (!spherical_image_map(x, y, z,Image, xcoor, ycoor)) 
  507.       return(1);
  508.     break;
  509.   case CYLINDRICAL_MAP:
  510.     if (!cylindrical_image_map(x, y, z,Image, xcoor, ycoor)) 
  511.       return(1);
  512.     break;
  513.   case TORUS_MAP:
  514.     if (!torus_image_map(x, y, z,Image, xcoor, ycoor)) 
  515.       return(1);
  516.     break;
  517.   default:
  518.     if (!planar_image_map(x, y, z,Image, xcoor, ycoor)) 
  519.       return(1);
  520.     break;
  521.   }
  522.  
  523.   /* Now make sure the point is on the image */
  524.   *ycoor += Small_Tolerance;
  525.   *xcoor += Small_Tolerance;
  526.   /* Compensate for y coordinates on the images being upsidedown */
  527.   *ycoor = (DBL)Image->iheight - *ycoor;
  528.  
  529.   if (*xcoor < 0.0)
  530.     *xcoor +=(DBL) Image->iwidth;
  531.   else if (*xcoor >=(DBL)Image->iwidth)
  532.     *xcoor -= (DBL)Image->iwidth;
  533.  
  534.   if (*ycoor < 0.0)
  535.     *ycoor += (DBL)Image->iheight;
  536.   else if (*ycoor >= (DBL)Image->iheight)
  537.     *ycoor -= (DBL)Image->iheight;
  538.  
  539.   if (Options & DEBUGGING)
  540.     printf ("\nmap %g %g %g xcoor %g ycoor %g ih %d iw %d\n", x, y, z, *xcoor, *ycoor,Image->iheight,Image->iwidth);
  541.  
  542.   if ((*xcoor >= (DBL)Image->iwidth) ||
  543.     (*ycoor >= (DBL)Image->iheight) ||
  544.     (*xcoor < 0.0) || (*ycoor < 0.0)) 
  545.     {
  546.     printf ("\nPicture index out of range\n");
  547.     close_all();
  548.     exit (1);
  549.     }
  550.   return(0);    
  551.   }
  552.  
  553. static void no_interpolation(Image, xcoor, ycoor, colour, index)
  554. IMAGE *Image;     
  555. DBL xcoor, ycoor;
  556. COLOUR *colour;
  557. int *index;
  558.   {  
  559.   struct Image_Line *line;
  560.   int iycoor, ixcoor;
  561.   IMAGE_COLOUR *map_colour;
  562.  
  563.   if (xcoor < 0.0)
  564.     xcoor += (DBL)Image->iwidth;
  565.   else if (xcoor >= (DBL)Image->iwidth)
  566.     xcoor -= (DBL)Image->iwidth;
  567.   if (ycoor < 0.0)
  568.     ycoor += (DBL)Image->iheight;
  569.   else if (ycoor >=(DBL) Image->iheight)
  570.     ycoor -= (DBL) Image->iheight;
  571.  
  572.   iycoor=(int)ycoor;
  573.   ixcoor=(int)xcoor;
  574.   if (Image->Colour_Map == NULL) 
  575.     {
  576.     line = &Image->data.rgb_lines[iycoor];
  577.     colour -> Red += (DBL) line->red[ixcoor]/255.0;
  578.     colour -> Green += (DBL) line->green[ixcoor]/255.0;
  579.     colour -> Blue += (DBL) line->blue[ixcoor]/255.0;
  580.     *index = -1;
  581.     } 
  582.   else 
  583.     {   
  584.     *index = Image->data.map_lines[iycoor][ixcoor];
  585.     map_colour = &Image->Colour_Map[*index];
  586.     /*printf ("icat index %d xc %d yc %d  CLR %d %d %d %d\n",*index,ixcoor,iycoor, 
  587.     map_colour->Red,map_colour->Green,map_colour->Blue,map_colour->Filter ); */
  588.     colour -> Red   += (DBL) map_colour->Red/255.0;
  589.     colour -> Green += (DBL) map_colour->Green/255.0;
  590.     colour -> Blue  += (DBL) map_colour->Blue/255.0;
  591.     colour -> Filter += (DBL) map_colour->Filter/255.0;
  592.     }
  593.  
  594.     if (Options & DEBUGGING)
  595.       printf ("\n no_interpolation index %d xc %d yc %d \n",*index,ixcoor,iycoor);
  596.   }
  597.  
  598. /* Interpolate color and filter values when mapping */
  599.  
  600. static void Interp(Image, xcoor, ycoor, colour, index)
  601. IMAGE *Image;     
  602. DBL xcoor, ycoor;
  603. COLOUR *colour;
  604. int *index;
  605.   {  
  606.   int iycoor, ixcoor,i;
  607.   int Corners_Index[4];
  608.   DBL Index_Crn[4];
  609.   COLOUR Corner_Colour[4];
  610.   DBL Red_Crn[4];
  611.   DBL Green_Crn[4];
  612.   DBL Blue_Crn[4];
  613.   DBL Filter_Crn[4];
  614.   DBL val1, val2, val3,val4;
  615.  
  616.   iycoor=(int)ycoor;
  617.   ixcoor=(int)xcoor;
  618.   for(i=0;i<4;i++)
  619.     {
  620.     Make_Colour(&Corner_Colour[i],0.0,0.0,0.0);
  621.     Corner_Colour[i].Filter=0.0; 
  622.     }
  623.   /* OK, now that you have the corners, what are you going to do with them? */
  624.   if(Image->Interpolation_Type==BILINEAR)
  625.     {
  626.     no_interpolation(Image,(DBL)ixcoor+1,(DBL)iycoor,&Corner_Colour[0],&Corners_Index[0]);
  627.     no_interpolation(Image,(DBL)ixcoor,(DBL)iycoor,&Corner_Colour[1],&Corners_Index[1]);
  628.     no_interpolation(Image,(DBL)ixcoor+1,(DBL)iycoor-1,&Corner_Colour[2],&Corners_Index[2]);
  629.     no_interpolation(Image,(DBL)ixcoor,(DBL)iycoor-1,&Corner_Colour[3],&Corners_Index[3]);
  630.     for(i=0;i<4;i++)
  631.       {
  632.       Red_Crn[i] = Corner_Colour[i].Red;
  633.       Green_Crn[i] = Corner_Colour[i].Green;
  634.       Blue_Crn[i] = Corner_Colour[i].Blue;
  635.       Filter_Crn[i] = Corner_Colour[i].Filter;
  636.       /* printf("Crn %d = %lf %lf %lf\n",i,Red_Crn[i],Blue_Crn[i],Green_Crn[i]); */
  637.       }
  638.  
  639.     val1 = bilinear(Red_Crn,xcoor,ycoor);
  640.     val2 = bilinear(Green_Crn,xcoor,ycoor);
  641.     val3 = bilinear(Blue_Crn,xcoor,ycoor);
  642.     val4 = bilinear(Filter_Crn,xcoor,ycoor);
  643.     }
  644.   if(Image->Interpolation_Type==NORMALIZED_DIST)
  645.     {
  646.     no_interpolation(Image,(DBL)ixcoor,(DBL)iycoor-1,&Corner_Colour[0],&Corners_Index[0]);
  647.     no_interpolation(Image,(DBL)ixcoor+1,(DBL)iycoor-1,&Corner_Colour[1],&Corners_Index[1]);
  648.     no_interpolation(Image,(DBL)ixcoor,(DBL)iycoor,&Corner_Colour[2],&Corners_Index[2]);
  649.     no_interpolation(Image,(DBL)ixcoor+1,(DBL)iycoor,&Corner_Colour[3],&Corners_Index[3]);
  650.     for(i=0;i<4;i++)
  651.       {
  652.       Red_Crn[i] = Corner_Colour[i].Red;
  653.       Green_Crn[i] = Corner_Colour[i].Green;
  654.       Blue_Crn[i] = Corner_Colour[i].Blue;
  655.       Filter_Crn[i] = Corner_Colour[i].Filter;
  656.       /* printf("Crn %d = %lf %lf %lf\n",i,Red_Crn[i],Blue_Crn[i],Green_Crn[i]); */
  657.       }
  658.  
  659.     val1 = norm_dist(Red_Crn,xcoor,ycoor);
  660.     val2 = norm_dist(Green_Crn,xcoor,ycoor);
  661.     val3 = norm_dist(Blue_Crn,xcoor,ycoor);
  662.     val4 = norm_dist(Filter_Crn,xcoor,ycoor);
  663.     }
  664.  
  665.   colour->Red   += val1;
  666.   colour->Green += val2;
  667.   colour->Blue  += val3;
  668.   colour->Filter += val4;
  669.   /* printf("Final = %lf %lf %lf\n",val1,val2,val3);  */
  670.   /* use bilinear for index try average later */
  671.   for(i=0;i<4;i++)
  672.     Index_Crn[i] = (DBL) Corners_Index[i];
  673.   if(Image->Interpolation_Type==BILINEAR)
  674.     {
  675.     *index = (int)(bilinear(Index_Crn,xcoor,ycoor)+0.5);
  676.     }
  677.   if(Image->Interpolation_Type==NORMALIZED_DIST)
  678.     {
  679.     *index = (int)(norm_dist(Index_Crn,xcoor,ycoor)+0.5);
  680.     }
  681.   }
  682.  
  683. /* These interpolation techniques are taken from an article by */
  684. /* Girish T. Hagan in the C Programmer's Journal V 9 No. 8 */
  685. /* They were adapted for POV-Ray by CdW */
  686.  
  687. static DBL bilinear (corners, x,y)
  688. DBL *corners, x,y;   
  689.   {
  690.   DBL p,q;
  691.   DBL val = 0.0;
  692.  
  693.   p = x - (int) x;
  694.   q = y - (int) y;
  695.   if ((p==0.0) && (q==0.0))
  696.     return(*corners); /* upper left */
  697.  
  698.   val = (p*q* *corners) + (q*(1-p)* *(corners+1)) + 
  699.   (p*(1-q)* *(corners+2)) + ((1-p)*(1-q)* *(corners+3));
  700.   return(val);
  701.   }  
  702.  
  703.  
  704. #define MAX_PTS 4
  705. #define PYTHAGOREAN_SQ(a,b)  ( (a)*(a) + (b)*(b) )
  706.  
  707. static DBL norm_dist(corners,x,y)
  708. DBL *corners,x,y;
  709.   {
  710.   register int i;
  711.  
  712.   DBL p,q;
  713.   DBL wts[MAX_PTS];
  714.   DBL sum_inv_wts = 0.0;
  715.   DBL sum_I  = 0.0;
  716.  
  717.   p = x - (int) x;
  718.   q = y - (int) y;
  719.  
  720.   if( (p==0.0) && (q==0.0) )
  721.     return(*corners); /* upper left */
  722.  
  723.   wts[0] = PYTHAGOREAN_SQ(p,q);
  724.   wts[1] = PYTHAGOREAN_SQ(1-p,q);
  725.   wts[2] = PYTHAGOREAN_SQ(p,1-q);
  726.   wts[3] = PYTHAGOREAN_SQ(1-p,1-q);
  727.  
  728.   for(i=0; i<MAX_PTS; i++)
  729.     {
  730.     sum_inv_wts += 1/wts[i];
  731.     sum_I += *(corners+i)/wts[i];
  732.     }
  733.  
  734.   return( sum_I /sum_inv_wts );
  735.   } 
  736.  
  737. IMAGE *Copy_Image (Old)
  738. IMAGE *Old;
  739.   {
  740.   return (Old);
  741.   }
  742.  
  743. IMAGE *Create_Image ()
  744.   {
  745.   IMAGE *Image;
  746.  
  747.   if ((Image = (IMAGE *)malloc(sizeof(IMAGE))) == NULL)
  748.     MAError("image file");
  749.   Image->File_Type = 0;
  750.   Image->Map_Type = PLANAR_MAP;
  751.   Image->Interpolation_Type = NO_INTERPOLATION;
  752.   Image->Once_Flag = FALSE;
  753.   Image->Use_Colour_Flag= FALSE;
  754.   Make_Vector (&Image->Gradient, 1.0, -1.0, 0.0);
  755.   return (Image);
  756.   }
  757.  
  758. void Destroy_Image (Image)
  759. IMAGE *Image;
  760.   {
  761.   }
  762.